/*-
 * Copyright (c) 2013-2014 Shang-Rong Cai. All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE AUTHORS AND CONTRIBUTORS ``AS IS'' AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED.  IN NO EVENT SHALL AUTHOR OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 * SUCH DAMAGE.
 *
 * $ sgmtest.c    Rev. 0.3.5    2014-06-19T16:24+08:00 $
 */


/*
 * 2D triangle grid in 3D coordinate of rectangle domain
 *
 * Domain: [x0, x1] by [y0, y1] (rectangle)
 * Partition into m parts in x-direction, n parts in y-direction
 * Grid: 2 * m * n elements
 */


/*
 * PDE:
 *     {-\triangle u = 0  in domain [-3.0, 5.0] by [-3.5, 2.5]
 *     { u(x, y) = 0.0    on boundary \Gamma _{0} = {x = -3.0, x = 5.0, y = 2.5}
 *     { u(x, y) = cos((x - 1.0) / 8.0 * pi)  on \Gamma _{1} = {y = -3.5}
 */


#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <stdint.h>
#include <errno.h>

#include <did.h>

#include <pregrid.h>
#include <trigrid.h>
#include <itermth.h>


/* Compute Stiffness Matrix */
int stiffness_matrix_triangle(int clz, int *vlst, double *coor, cidk_t *cidk);

int local_stiffness_triangle(int clz, int *vlst, double *coor);


int find_position_recursive(int trg, int lstz, cidbl_t *ilst, cidbl_t **addr, int *flg);

int test_find_index(int ridx, int cidx, cidk_t *cidk);

int cidk_value_add(int ridx, int cidx, double value, cidk_t *cidk);


/* Save DID File */
int save_did_i32_f64(FILE *fpof, uint64_t fsz, int *fst_idx, uint64_t ssz, int *snd_idx, uint64_t esz, double *dat_ent);

int save_matrix_did(const char *filename, cidk_t *cidk);


/* Boundary Condition Nodeset */
typedef struct
{
    size_t ndz;
    int *idx_lst;
} bcnd_t;

int bcnd_init(size_t nz, bcnd_t *bc);

int bcnd_dest(bcnd_t *bc);

int bcnd_stdout(const bcnd_t *bc);


/* Dirichlet Boundary Condition For 0 (Others) */
int domain_dirichlet_bc_0_create(trigrid_t *tgrd, double x0, double x1, double y0, double y1, double epsilon, bcnd_t *bcnd);

#define domain_dirichlet_bc_0_delete bcnd_dest

/* Dirichlet Boundary Condition For  cos((x - 1.0) / 8.0 * pi) (Others) */
/*
 * (x0 + x1) / 2.0 == 1.0
 * (x1 - x0) == 8.0
 */
int domain_dirichlet_bc_1_create(trigrid_t *tgrd, double x0, double x1, double y0, double y1, double epsilon, bcnd_t *bcnd);

#define domain_dirichlet_bc_1_delete bcnd_dest


/* Boundary Condition Matrix */
int set_boundary_dirichlet_matrix(const bcnd_t *bcnd, cidk_t *cidk);


int vector_save_vtk_append(const char *filename, const size_t len, double *vec);


/* Boundary Condition For Right Hand Side Vector */
int boundary_0_function(const trigrid_t *tgrd, const bcnd_t *bcnd, double *vec);

int boundary_1_function(const trigrid_t *tgrd, const bcnd_t *bcnd, double *vec);


int main(int argc, char *argv[])
{
    int ierr = ERRNO_UNDEFINED;

    trigrid_t tgrd;  /* Triangular mesh in rectangular domain */

    int rctr = 2;  /* Square divide as two triangle */
    int trvz = 3;  /* Each triangle has three vertices */

    double x0, x1;  /* Rectangle region of x-coordinate start and end */
    double y0, y1;  /* Rectangle region of y-coordinate start and end */

    int m, n;  /* Partitions for x and y directions */

    cidk_t cidk;  /* Stiffness matrix */

    double epsilon;  /* Small Shift For Boundary Conditions */

    bcnd_t bc_d_0;  /* Dirichlet boundary conditions 0 */
    bcnd_t bc_d_1;  /* Dirichlet boundary conditions 1 */

    double *vec_rhs;  /* Right hand side of Ax=b */
    double *vec_sol;  /* Solution of Ax=b */

    int buf_iter;  /* Iterations of method */

    /* Set domain */
    x0 = -3.0;
    x1 = 5.0;

    y0 = -3.5;
    y1 = 2.5;

    m = 8;
    n = 6;

    /* Generate Vertices Coordinates And Cell */
    ierr = trigrid_create(x0, x1, y0, y1, m, n, &tgrd);

    /* STDOUT tgrd */
#ifdef DEBUG_INFO
    ierr = trigrid_stdout(&tgrd);
#endif

    /* Save Doamin Geometry To A VTK File */
    ierr = trigrid_save_vtk(&tgrd, "sample.vtk");

    /* Create Stiffness Matrix */
    ierr = trigrid_create_cidk(&tgrd, &cidk);

    /* Stiffness Matrix */
    ierr = stiffness_matrix_triangle(tgrd.elmz, tgrd.vlst, tgrd.coor, &cidk);
#ifdef DEBUG_INFO
    ierr = cidk_stdout(&cidk);
#endif

    /* Save Stiffness Matrix (Without Setting Boundary Condition) */
    ierr = save_matrix_did("sample.did", &cidk);

    /* Create Boundary Node Set */
    epsilon = (x1 - x0) / m / 2;
#ifdef DEBUG_INFO
    printf("epsilon = %lE\n", epsilon);
#endif

    ierr = domain_dirichlet_bc_0_create(&tgrd, x0, x1, y0, y1, epsilon, &bc_d_0);
    ierr = domain_dirichlet_bc_1_create(&tgrd, x0, x1, y0, y1, epsilon, &bc_d_1);

#ifdef DEBUG_INFO
    ierr = bcnd_stdout(&bc_d_0);
    ierr = bcnd_stdout(&bc_d_1);
#endif

    /* Set Boundary Condition To Matrix */
    ierr = set_boundary_dirichlet_matrix(&bc_d_0, &cidk);
    ierr = set_boundary_dirichlet_matrix(&bc_d_1, &cidk);

    /* Save Stiffness Matrix (Set Boundary Condition) */
    ierr = save_matrix_did("sample_mtx.did", &cidk);

    /* Create Right Hand Side And Solution Vectors */
    vec_rhs = (double*) calloc(tgrd.vtxz, sizeof(double));
    vec_sol = (double*) calloc(tgrd.vtxz, sizeof(double));

    ierr = vector_set_zeros(tgrd.vtxz, vec_rhs);
    ierr = vector_set_zeros(tgrd.vtxz, vec_sol);

    /* Set Boundary Condition To Right Hand Side */
    ierr = boundary_1_function(&tgrd, &bc_d_1, vec_rhs);
    ierr = boundary_0_function(&tgrd, &bc_d_0, vec_rhs);
#ifdef DEBUG_INFO
    ierr = vector_stdout(tgrd.vtxz, vec_rhs);
#endif

    /* Iterative Method Solve Ax=b And Save Solution */
    printf("\nJacobi Method\n");
    ierr = jacobi_cidk(&cidk, vec_rhs, 1000, 1.0e-10, 1.0e-5, 1.0, vec_sol, &buf_iter);
    ierr = vector_save_vtk_append("sample_solution.jac10", tgrd.vtxz, vec_sol);

    printf("\nWeighted Jacobi Method (Omega = 0.6)\n");
    ierr = vector_set_zeros(tgrd.vtxz, vec_sol);
    ierr = jacobi_cidk(&cidk, vec_rhs, 1000, 1.0e-10, 1.0e-5, 0.8, vec_sol, &buf_iter);
    ierr = vector_save_vtk_append("sample_solution.jac06", tgrd.vtxz, vec_sol);

    printf("\nGauss Seidel Method\n");
    ierr = vector_set_zeros(tgrd.vtxz, vec_sol);
    ierr = sor_cidk(&cidk, vec_rhs, 1000, 1.0e-10, 1.0e-5, 1.0, vec_sol, &buf_iter);
    ierr = vector_save_vtk_append("sample_solution.sor10", tgrd.vtxz, vec_sol);

    printf("\nSOR Method (Omega = 0.5)\n");
    ierr = vector_set_zeros(tgrd.vtxz, vec_sol);
    ierr = sor_cidk(&cidk, vec_rhs, 1000, 1.0e-10, 1.0e-5, 0.5, vec_sol, &buf_iter);
    ierr = vector_save_vtk_append("sample_solution.sor05", tgrd.vtxz, vec_sol);

    printf("\nSOR Method (Omega = 1.5)\n");
    ierr = vector_set_zeros(tgrd.vtxz, vec_sol);
    ierr = sor_cidk(&cidk, vec_rhs, 1000, 1.0e-10, 1.0e-5, 1.5, vec_sol, &buf_iter);
    ierr = vector_save_vtk_append("sample_solution.sor15", tgrd.vtxz, vec_sol);

    /* Free */
    ierr = trigrid_delete(&tgrd);

    ierr = trigrid_delete_cidk(&cidk);

    ierr = domain_dirichlet_bc_0_delete(&bc_d_0);
    ierr = domain_dirichlet_bc_0_delete(&bc_d_1);

    free(vec_rhs);
    free(vec_sol);

    return ierr;
}


int stiffness_matrix_triangle(int clz, int *vlst, double *coor, cidk_t *cidk)
{
    int ierr = ERRNO_UNDEFINED;
    int ii, ij, ik;

    refelm_t refelm;

    int vidx[3];  /* Vertex index (physical element) */

    double p1;  /* x_1 - x_0 */
    double p2;  /* x_2 - x_0 */
    double q1;  /* y_1 - y_0 */
    double q2;  /* y_2 - y_0 */

    double det_ma;  /* Determinate: p1 * q2 - q1 * p2; (absolute value is area of triangle element) */

    double emtx[9];  /* Local stiffness matrix */

    double dpsi_dx[3], dpsi_dy[3]; /* Partial derivative for each basis (physical) */

    /* Test Find Index (Part of Computing Stiffness Matrix) */
#ifdef DEBUG_INFO
    ierr = test_find_index(1, 1, cidk);
    ierr = test_find_index(10, 14, cidk);
    ierr = test_find_index(84, 83, cidk);
#endif

    /* Reference Element */
    ierr = refelm_triangle_init(&refelm);

    /* Local Stiffness Matrix And Assembling */
    for (ik = 0; ik < clz; ik++)
    {
        vidx[0] = vlst[ik * 3];
        vidx[1] = vlst[ik * 3 + 1];
        vidx[2] = vlst[ik * 3 + 2];

    #ifdef DEBUG_INFO
        printf("%3d: ", ik);
        printf("%3d: (% .3le, % .3le),\n", vidx[0], coor[vidx[0] * 3], coor[vidx[0] * 3 + 1]);
        printf("     %3d: (% .3le, % .3le),\n", vidx[1], coor[vidx[1] * 3], coor[vidx[1] * 3 + 1]);
        printf("     %3d: (% .3le, % .3le).\n", vidx[2], coor[vidx[2] * 3], coor[vidx[2] * 3 + 1]);
    #endif

        p1 = coor[vidx[1] * 3] - coor[vidx[0] * 3];
        p2 = coor[vidx[2] * 3] - coor[vidx[0] * 3];

        q1 = coor[vidx[1] * 3 + 1] - coor[vidx[0] * 3 + 1];
        q2 = coor[vidx[2] * 3 + 1] - coor[vidx[0] * 3 + 1];

        det_ma = p1 * q2 - q1 * p2;

    #ifdef DEBUG_INFO
        printf("     det(A) = % 7le\n", det_ma);
    #endif
        /* Partial Derivative For Each Basis (Physical) */
        dpsi_dx[0] = q2 * refelm.refpvr[0] - q1 * refelm.refpvs[0];
        dpsi_dy[0] = p1 * refelm.refpvs[0] - p2 * refelm.refpvr[0];

        dpsi_dx[1] = q2 * refelm.refpvr[1] - q1 * refelm.refpvs[1];
        dpsi_dy[1] = p1 * refelm.refpvs[1] - p2 * refelm.refpvr[1];

        dpsi_dx[2] = q2 * refelm.refpvr[2] - q1 * refelm.refpvs[2];
        dpsi_dy[2] = p1 * refelm.refpvs[2] - p2 * refelm.refpvr[2];

        for (ii = 0; ii < 3; ii++)
        {
        #ifdef DEBUG_INFO
            printf("[");
        #endif
            for (ij = 0; ij < 3; ij++)
            {
                emtx[ii * 3 + ij] = (dpsi_dx[ii] * dpsi_dx[ij] + dpsi_dy[ii] * dpsi_dy[ij]) / fabs(det_ma) * refelm.qdtwgt[0];
            #ifdef DEBUG_INFO
                printf(" % .7lE", emtx[ii * 3 + ij]);
            #endif
                ierr = cidk_value_add(vidx[ii], vidx[ij], emtx[ii * 3 + ij], cidk);
            }
        #ifdef DEBUG_INFO
            printf("]\n");
        #endif
        }
    }

    /* Finalize */
    ierr = refelm_triangle_dest(&refelm);

    return ierr;
}


int find_position_recursive(int trg, int lstz, cidbl_t *ilst, cidbl_t **addr, int *flg)
{
    int ierr = ERRNO_UNDEFINED;
    int half;

    *flg = 0;

    if (lstz <= 0)
    {
    #ifdef DEBUG_INFO
        printf("Branch 0:\n");
    #endif
        perror("Length should not be zero!\n");
        ierr = ERRNO_OUT_OF_RANGE_SIZE;
    }
    else if (lstz == 1)
    {
    #ifdef DEBUG_INFO
        printf("Branch 1:\n");
    #endif
        if (trg == ilst->index)
        {
            *addr = ilst;
            *flg = 1;
        }
        else
        {
            printf("The ilst->index is %d.\n", ilst->index);
            printf("Index %d not found!\n", trg);
        }
    }
    else if (lstz == 2)
    {
    #ifdef DEBUG_INFO
        printf("Branch 2:\n");
    #endif
        if (trg == ilst[0].index)
        {
            *addr = ilst;
            *flg = 1;
        }
        else if (trg == ilst[1].index)
        {
            *addr = &ilst[1];
            *flg = 1;
        }
        else
        {
            printf("Index %d not found!\n", trg);
        }
    }
    else
    {
    #ifdef DEBUG_INFO
        printf("Branch 3:\n");
    #endif
        half = lstz / 2;
        if (trg < ilst[half].index)
        {
            ierr = find_position_recursive(trg, half, ilst, addr, flg);
        }
        else if (trg == ilst[half].index)
        {
            *addr = &ilst[half];
            *flg = 1;
        }
        else if (trg > ilst[half].index)
        {
            if (lstz % 2 == 0)
            {
                ierr = find_position_recursive(trg, half, &ilst[half], addr, flg);
            }
            else
            {
                ierr = find_position_recursive(trg, half, &ilst[half + 1], addr, flg);
            }
        }
    }

    EXIT_FIND_POSITION:

    return ierr;
}


int test_find_index(int ridx, int cidx, cidk_t *cidk)
{
    int ierr = ERRNO_UNDEFINED;
    int lstz;
    cidbl_t *addr;
    int flag;

    lstz = cidk->sptr[ridx].tail - cidk->sptr[ridx].head;

    ierr = find_position_recursive(cidx, lstz, &cidk->item[cidk->sptr[ridx].head], &addr, &flag);

    if (flag == 1)
    {
        printf("Found it! (%d, %d): % .7lE\n", ridx, cidx, addr->value);
    }

    return ierr;
}


int cidk_value_add(int ridx, int cidx, double value, cidk_t *cidk)
{
    int ierr = ERRNO_UNDEFINED;
    int lstz;
    cidbl_t *addr;
    int flag;

    lstz = cidk->sptr[ridx].tail - cidk->sptr[ridx].head;

    ierr = find_position_recursive(cidx, lstz, &cidk->item[cidk->sptr[ridx].head], &addr, &flag);

    if (flag == 1)
    {
        addr->value += value;
    }
    else
    {
        printf("Position (%d, %d) not found!\n", ridx, cidx);
    }

    return ierr;
}


int bcnd_init(size_t nz, bcnd_t *bc)
{
    int ierr = ERRNO_UNDEFINED;

    bc->ndz = nz;

    bc->idx_lst = calloc(bc->ndz, sizeof(int));

    return ierr;
}

int bcnd_dest(bcnd_t *bc)
{
    int ierr = ERRNO_UNDEFINED;

    free(bc->idx_lst);

    return ierr;
}


int bcnd_stdout(const bcnd_t *bc)
{
    int ierr = ERRNO_UNDEFINED;

    size_t zi;

    int flg_nl;  /* Flag new line; 0: no newline, 1: newline printed */

    printf("Boundary Condition Node Set:");

    for (zi = 0; zi < bc->ndz; zi++)
    {
        if (zi % 20 == 0)
        {
            printf("\n");
        }
        printf(" %3d", bc->idx_lst[zi]);

        if (zi / 20 == 19)
        {
            flg_nl = 1;
        }
        else
        {
            flg_nl = 0;
        }
    }

    if (flg_nl == 0)
    {
        printf("\n");
    }

    return ierr;
}


/* Dirichlet Boundary Condition For 0 (Others) */
int domain_dirichlet_bc_0_create(trigrid_t *tgrd, double x0, double x1, double y0, double y1, double epsilon, bcnd_t *bcnd)
{
    int ierr = ERRNO_UNDEFINED;

    size_t zi;
    size_t tb_idx;  /* Test index of vertex coordinate block */
    size_t buf_ndz;  /* Number of nodes */

    double dx_l;  /* Distance to x0 */
    double dx_r;  /* Distance to x1 */
    double dy_b;  /* Distance to y0 */
    double dy_t;  /* Distance to y1 */

    double cur_cx, cur_cy;

    buf_ndz = 0;

    for (zi = 0; zi < tgrd->vtxz; zi++)
    {
        tb_idx = zi * tgrd->sdim;

        cur_cx = tgrd->coor[tb_idx];
        cur_cy = tgrd->coor[tb_idx + 1];

        dx_l = fabs(cur_cx - x0);
        dx_r = fabs(cur_cx - x1);

        dy_b = fabs(cur_cy - y0);  /* For other boundary condition */
        dy_t = fabs(cur_cy - y1);

        if (dx_l < epsilon || dx_r < epsilon || dy_t < epsilon)
        {
            buf_ndz++;
        #ifdef DEBUG_INFO
            printf("%4zd: Index: %4zd, %lf, %lf, %lf, %lf\n", buf_ndz, zi, dx_l, dx_r, dy_b, dy_t);
        #endif
        }
    }

    /* Create Boundary Node Set */
    ierr = bcnd_init(buf_ndz, bcnd);

    buf_ndz = 0;

    for (zi = 0; zi < tgrd->vtxz; zi++)
    {
        tb_idx = zi * tgrd->sdim;

        cur_cx = tgrd->coor[tb_idx];
        cur_cy = tgrd->coor[tb_idx + 1];

        dx_l = fabs(cur_cx - x0);
        dx_r = fabs(cur_cx - x1);

        dy_b = fabs(cur_cy - y0);  /* For other boundary condition */
        dy_t = fabs(cur_cy - y1);

        if (dx_l < epsilon || dx_r < epsilon || dy_t < epsilon)
        {
            bcnd->idx_lst[buf_ndz] = (int) zi;
            buf_ndz++;
        }
    }

    return ierr;
}


/* Dirichlet Boundary Condition For  cos((x - 1.0) / 8.0 * pi) (Others) */
/*
 * (x0 + x1) / 2.0 == 1.0
 * (x1 - x0) == 8.0
 */
int domain_dirichlet_bc_1_create(trigrid_t *tgrd, double x0, double x1, double y0, double y1, double epsilon, bcnd_t *bcnd)
{
    int ierr = ERRNO_UNDEFINED;

    size_t zi;
    size_t tb_idx;  /* Test index of vertex coordinate block */
    size_t buf_ndz;  /* Number of nodes */

    double dx_l;  /* Distance to x0 */
    double dx_r;  /* Distance to x1 */
    double dy_b;  /* Distance to y0 */
    double dy_t;  /* Distance to y1 */

    double cur_cx, cur_cy;

    buf_ndz = 0;

    for (zi = 0; zi < tgrd->vtxz; zi++)
    {
        tb_idx = zi * tgrd->sdim;

        cur_cx = tgrd->coor[tb_idx];
        cur_cy = tgrd->coor[tb_idx + 1];

        dx_l = fabs(cur_cx - x0);  /* For other boundary condition */
        dx_r = fabs(cur_cx - x1);  /* For other boundary condition */

        dy_b = fabs(cur_cy - y0);
        dy_t = fabs(cur_cy - y1);  /* For other boundary condition */

        if (dy_b < epsilon)
        {
            buf_ndz++;
        #ifdef DEBUG_INFO
            printf("%4zd: Index: %4zd, %lf, %lf, %lf, %lf\n", buf_ndz, zi, dx_l, dx_r, dy_b, dy_t);
        #endif
        }
    }

    /* Create Boundary Node Set */
    ierr = bcnd_init(buf_ndz, bcnd);

    buf_ndz = 0;

    for (zi = 0; zi < tgrd->vtxz; zi++)
    {
        tb_idx = zi * tgrd->sdim;

        cur_cx = tgrd->coor[tb_idx];
        cur_cy = tgrd->coor[tb_idx + 1];

        dx_l = fabs(cur_cx - x0);  /* For other boundary condition */
        dx_r = fabs(cur_cx - x1);  /* For other boundary condition */

        dy_b = fabs(cur_cy - y0);
        dy_t = fabs(cur_cy - y1);  /* For other boundary condition */

        if (dy_b < epsilon)
        {
            bcnd->idx_lst[buf_ndz] = (int) zi;
            buf_ndz++;
        }
    }

    return ierr;
}


int set_boundary_dirichlet_matrix(const bcnd_t *bcnd, cidk_t *cidk)
{
    int ierr = ERRNO_UNDEFINED;

    size_t zi;
    sidx_t si;

    sidx_t trgt_idx;

    for (zi = 0; zi < bcnd->ndz; zi++)
    {
        trgt_idx = bcnd->idx_lst[zi];
    #ifdef DEBUG_INFO
        printf("Target Index: %4d (Row)\n", trgt_idx);
    #endif

        for (si = cidk->sptr[trgt_idx].head; si < cidk->sptr[trgt_idx].tail; si++)
        {
        #ifdef DEBUG_INFO
            printf("    Column: %4d, Value: % le\n", cidk->item[si].index, cidk->item[si].value);
        #endif

            if (trgt_idx == cidk->item[si].index)
            {
                cidk->item[si].value = 1.0;
            }
            else
            {
                cidk->item[si].value = 0;
            }
        #ifdef DEBUG_INFO
            printf("    Column: %4d, Value: % le\n", cidk->item[si].index, cidk->item[si].value);
        #endif
        }
    }

    return ierr;
}


int vector_save_vtk_append(const char *filename, const size_t len, double *vec)
{
    int ierr = ERRNO_UNDEFINED;
    size_t zi;

    FILE *fpof;

    fpof = fopen(filename, "w");

    fprintf(fpof, "POINT_DATA %zd\n", len);
    fprintf(fpof, "SCALARS solution float 1\n");
    fprintf(fpof, "LOOKUP_TABLE default\n");

    for (zi = 0; zi < len; zi++)
    {
        fprintf(fpof, " % .7le\n", vec[zi]);
    }

    fclose(fpof);

    return ierr;
}


/*
 * u(x, y) = 0.0  on boundary \Gamma _{0} = {x = -3.0, x = 5.0, y = 2.5}
 */
int boundary_0_function(const trigrid_t *tgrd, const bcnd_t *bcnd, double *vec)
{
    int ierr = ERRNO_UNDEFINED;
    size_t zi;

    for (zi = 0; zi < bcnd->ndz; zi++)
    {
        vec[bcnd->idx_lst[zi]] = 0;
    }

    return ierr;
}


/*
 * { u(x, y) = cos((x - 1.0) / 8.0 * pi)  on \Gamma _{1} = {y = -3.5}
 */
int boundary_1_function(const trigrid_t *tgrd, const bcnd_t *bcnd, double *vec)
{
    int ierr = ERRNO_UNDEFINED;

    size_t zi;
    size_t nd_bki;  /* Node index of vertex coordinate block */

    double xv;  /* x-coordinate */
    double bv;  /* Boundary value */

    const double local_pi = 3.14159265358;

    for (zi = 0; zi < bcnd->ndz; zi++)
    {
        nd_bki = (size_t) bcnd->idx_lst[zi] * tgrd->sdim;

        xv = tgrd->coor[nd_bki];
        bv = cos((xv - 1.0) / 8.0 * local_pi);

    #ifdef DEBUG_INFO
        printf("Boundary Condition Node Index: %4d; x-coordinate % le\n", bcnd->idx_lst[zi], xv);
        printf("Boundary Value: % le\n", bv);
    #endif

        vec[bcnd->idx_lst[zi]] = bv;
    }

    return ierr;
}


int save_did_i32_f64(FILE *fpof, uint64_t fsz, int *fst_idx, uint64_t ssz, int *snd_idx,
uint64_t esz, double *dat_ent)
{
    int ierr = ERRNO_UNDEFINED;

    did_head_t hdprt;
    did_meta_t mtprt;
    did_body_t bdprt;

    ierr = did_head_set_info(&hdprt, DID_HEAD_INDEX_T_I32, DID_HEAD_DATA_T_F64);
    ierr = did_meta_set_info(&mtprt, &hdprt, fsz, ssz, esz);

    ierr = did_body_create(&bdprt, &hdprt, &mtprt, (void*) fst_idx, (void*) snd_idx, (void*) dat_ent);

    ierr = did_write(fpof, &hdprt, &mtprt, &bdprt);

    ierr = did_body_delete(&bdprt);

    return ierr;
}


int save_matrix_did(const char *filename, cidk_t *cidk)
{
    int ierr = ERRNO_UNDEFINED;

    size_t zi;
    sidx_t si;

    size_t cmpidxz;  /* Compressed index size */
    size_t cnt;

    int32_t *fst_idx;
    int32_t *snd_idx;
    double *dat_ent;

    FILE *fpof;

    /* Check Scheme */
    switch (cidk->flag)
    {
        case ROW_WISE:
            cmpidxz = cidk->r_sz;
            break;
        case COL_WISE:
            cmpidxz = cidk->c_sz;
            break;
        case RDR_NONE:
            ierr = ERRNO_SUPPORT_NONE;
            goto EXIT_SAVE_MATRIX_DID;
        default:
            ierr = ERRNO_UNKNOWN_TYPE;
            goto EXIT_SAVE_MATRIX_DID;
            break;
    }

    /* Allocate Memory */
    fst_idx = (int32_t*) calloc(cidk->e_sz, sizeof(int32_t));
    snd_idx = (int32_t*) calloc(cidk->e_sz, sizeof(int32_t));
    dat_ent = (double*) calloc(cidk->e_sz, sizeof(double));

    cnt = 0;

    if (cidk->flag == ROW_WISE)
    {
        for (zi = 0; zi < cmpidxz; zi++)
        {
            for (si = cidk->sptr[zi].head; si < cidk->sptr[zi].tail; si++)
            {
                fst_idx[cnt] = (int32_t) zi;
                snd_idx[cnt] = (int32_t) cidk->item[si].index;
                dat_ent[cnt] = cidk->item[si].value;
                cnt++;
            }
        }
    }
    else
    {
        for (zi = 0; zi < cmpidxz; zi++)
        {
            for (si = cidk->sptr[zi].head; si < cidk->sptr[zi].tail; si++)
            {
                fst_idx[cnt] = (int32_t) cidk->item[si].index;
                snd_idx[cnt] = (int32_t) zi;
                dat_ent[cnt] = cidk->item[si].value;
                cnt++;
            }
        }
    }

    /* Write To File */
    fpof = fopen(filename, "w");
    if (fpof == NULL)
    {
        perror(filename);
        ierr = errno;
        goto EXIT_SAVE_MATRIX_DID_FREE_MEMORY;
    }

    ierr = save_did_i32_f64(fpof, (uint64_t) cidk->r_sz, fst_idx, (uint64_t) cidk->c_sz, snd_idx, (uint64_t) cidk->e_sz, dat_ent);

    fclose(fpof);

    EXIT_SAVE_MATRIX_DID_FREE_MEMORY:
    free(fst_idx);
    free(snd_idx);
    free(dat_ent);

    EXIT_SAVE_MATRIX_DID:

    return ierr;
}
